Nonlinear steepest descent and the numerical solution of 
Riemann-Hilbert problems 

Sheehan Olver 1 and Thomas Trogdon 2 

1 School of Mathematics and Statistics 
The University of Sydney 
NSW 2006, Australia 

2 Department of Applied Mathematics 
University of Washington 

Campus Box 352420 
Seattle, WA, 98195, USA 

March 4, 2013 



Abstract 

The effective and efficient numerical solution of Riemann-Hilbert problems has been demonstrated 
in recent work. With the aid of ideas from the method of nonlinear steepest descent for Riemann- 
Hilbert problems, the resulting numerical methods have been shown numerically to retain accuracy as 
values of certain parameters become arbitrarily large. Remarkably, this numerical approach does not 
require knowledge of local parametrices; rather, the deformed contour is scaled near stationary points 
at a specific rate. The primary aim of this paper is to prove that this observed asymptotic accuracy is 
indeed achieved. To do so, we first construct a general theoretical framework for the numerical solution of 
Riemann-Hilbert problems. Second, we demonstrate the precise link between nonlinear steepest descent 
and the success of numerics in asymptotic regimes. In particular, we prove sufficient conditions for 
numerical methods to retain accuracy. Finally, we compute solutions to the homogeneous Painleve II 
equation and the modified Korteweg-de Vries equations to explicitly demonstrate the practical validity 
of the theory. 

1 Introduction 

Matrix-valued Riemann-Hilbert problems (RHPs) are of profound importance in modern applied analysis. 
In inverse scattering theory, solutions to the nonlinear Schrodinger equation, Korteweg de- Vries equation 
(KdV), Kadomtsev-Petviashvili I equation and many other integrable solutions can be written in terms of 
solutions of RHPs pQ. Orthogonal polynomials can also be rewritten in terms of solutions of RHPs [3]. 
The asymptotics of orthogonal polynomials is crucial for determining the distribution of eigenvalues of large 
random matrix ensembles [3]. In each of these applications, RHPs fulfill the role that integral representations 
play in classical asymptotic analysis. 

The way in which these RHPs are analyzed is through the method of nonlinear steepest descent [I]. 
RHPs can be deformed in the complex plane in much the same way as contour integrals. This allows the 
oscillatory nature of the problem to be changed to exponential decay. These deformed problems, which 
depend on a parameter, are solved approximately. The approximate solution is found explicitly and the 
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difference between the actual solution and this approximate solution tends to zero as a parameter becomes 
large. 

The method of nonlinear steepest descent was adapted by the authors in |15j for numerical purposes. 
We developed a method to reliably solve the Cauchy initial-value problem for KdV and modified KdV 
for all values of x and t. A benefit of the approach was that, unlike the standard method of nonlinear 
steepest descent, we did not require the knowledge of difficult-to-derive local parametrices. Instead, an 
approach was developed based on scaling contours. Convergence of the approximation was demonstrated 
through numerical tests, which further showed that the accuracy of the approximation actually increased 
in the asymptotic regime. The focus of the current paper is to derive sufficient conditions (which are often 
satisfied) for which we can prove that the approach of solving RHPs numerically on scaled contours will be 
guaranteed to be accurate in asymptotic regimes. We refer to this type of behavior as asymptotic stability 
or uniform approximation. 

In addition, we show the deep connection between the success of the numerical method and the success of 
the method of nonlinear steepest descent [4] . A notable conclusion is that one can expect that whenever the 
method of nonlinear steepest descent produces an asymptotic formula, the numerical method can be made 
asymptotically stable. Achieving this requires varying amounts of preconditioning of the RHP. This can vary 
from not deforming the RHP at all, all the way to using the full deformation needed by the analytical method. 
An important question is: "when can we stop deformations and a have a reliable numerical method?" Our 
main results are in <|5]and allow us to answer this question. In short, although we do not require the knowledge 
of local parametrices to construct the numerical method, their existence ensures that the numerical method 
remains accurate, and their explicit knowledge allows us to analyze the error of the approximation directly. 
In the last two sections we provide useful examples of where this arises in applications. 

The paper is structured as follows. We begin with some background material, followed by the precise 
definition of a RHP along with properties of an associated singular integral operator. This allows us, amongst 
other things, to address the regularity of solutions. Then, we use an abstract framework for the numerical 
solution of RHP which will allow us to address asymptotic accuracy in a more concise way. Additionally, 
other numerical methods (besides the one used for the applications) may fit within the framework of this 
paper. We review the philosophy and analysis behind nonlinear steepest descent and how it relates to 
our numerical framework. Then, we prove our main results which provide sufficient condition for uniform 
approximation. The numerical approach of is placed within the general framework, along with necessary 
assumptions which allows a realization of uniform approximation. We apply the theory to two RHPs. The 
first is a RHP representation of the homogeneous Painleve II transcendent 

u xx — xu — 2u 3 = 0, (1-1) 

for specific initial conditions and x < 0. The second is a RHP representation of the modified Korteweg-de 
Vries equation 

u t - 6u 2 u x + u xxx = 0, (1.2) 
for smooth exponentially decaying initial data in the so-called Painleve region |S]. 

2 Background Material 

We use this section to fix notation that will be used throughout the remainder of the manuscript. We reserve 
C and Ci to denote generic constants, and X and Y for Banach spaces. We denote the norm on a space X 
by || ■ || x- The notation C(X, Y) is used to denote the Banach algebra of all bounded linear operators from 
X to Y , When X = Y we write C(X) to simplify notation. The following lemma is of great use 

Lemma 2.1. Assume L 6 £(A, Y) has a bounded inverse L^ 1 G C(Y,X). Assume M £ C(X,Y) satisfies 

l|M " L|l< ^' 
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Then M is invertible ana 



Furthe 



IM-'W < n ^rr [7 • (2.1) 



||X^-M-i< "V 112 ' 1 ^^' (2.2) 
11 11 - 1 - IIL-^HIL - M|| V ; 

In what follows we are interested in functions denned on oriented contours in C. Assume T is piece- 
wise smooth, oriented and (most likely) self-intersecting. 70 is used to denote the set of self-intersections. 
Decompose r = Ti U • • • U Ti to its smooth, non-self-intersecting components. Define C^(Ti) to be the 
space of infinitely differentiable functions with compact support in I\, and C(Ti) to be the Banach space of 
continuous functions with the uniform norm. Define the space 



L 2 {T) = j/ measurable : ^ J |/(fc)| 2 |dfc| < 00 j 



We define, D, the distributional differentiation operator for functions defined on T \ 70. For a function 
<p £ C^°(Ti) we represent a linear function g via the dual pairing 

9{<P) = (sS^ry 

To Dig we associate the functional 

(ff,^)rv 

For / £ L 2 (T) consider = f\^ i , the restriction of / to Fj. In the case that the distribution Difi corresponds 
to a locally integrable function 

{Dif u <p) Tt = J Difi(kMk)dk = J fi(k)ip'(k)dk, 

for each i, we define 

Df{k) = Difi(k) iffcGr i \ 7o . 

This allows us to define 

H k (T) = {f£ L 2 (T) : Dif £ L 2 (T), j = 0,...,k}. 

We write W k, °°(T) for the Sobolev space with the L? norm replaced with the L°° norm. An important note 
is that we will be dealing with matrix-valued functions, and hence the definitions of all these spaces must 
be suitably extended. Since all finite-dimensional norms are equivalent, we can use the above definitions 
in conjunction with any matrix norm to define a norm for matrix-valued functions provided the norm is 
sub-additive. 



3 Theory of Riemann— Hilbert Problems in L 2 (T) 

Loosely speaking, a Riemann-Hilbert problem (RHP) is a boundary-value problem in the complex plane: 

Problem 3.1. Given an oriented contour T C C and a jump matrix G : T — > C 2x2 , find a bounded 
function $:C\r— >C 2x2 which is analytic everywhere in the complex plane except on T, such that 

<& + (z)=$-(z)G(z), forz£T,and (3.1) 
$(00) = I, (3.2) 

where $ + denotes the limit of $ as z approaches T from the left, $~ denotes the limit of $ as z approaches 
r from the right, and $(00) = lim| z |_j. 00 <&{z). We denote this RHP by [G; T]. 
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The definition is not sufficiently precise to compute solutions. We say that 4> is an L 2 solution of [G; T] 
normalized at infinity, if 

$(z)=/+~ / U ^ds = I + C r u{z), 

ZlTl J T S — Z 

for some u € L 2 (T) and $ + (z) = Q~ (z)G(z), for almost every zeT. 

Definition 3.1. The Cauchy transform pair is the pair of boundary values of the Cauchy integral Cr' 

C±u(z) = (C r w(z)) ± . (3.3) 

It is well known that for fairly general contours (including every contour we considered here) that the 
limits in (3.3) exist non-tangentially almost everywhere. The operators Cp are bounded from L 2 (T) to itself 
and satisfy the operator identity [5] 

C+ - Cf = I. (3.4) 

Remark 3.1. When deriving a RHP one must show that $ — / 6 ranCr- This is generally done using Hardy 
space methods /<§, 61. Brevity is the governing motivation for the above definition of a RHP. 

We convert the RHP into an equivalent singular integral equation (SIE). Assume $(z) = I + Cru(z) and 
substitute into (3.1) to obtain 

I + C£u(z) =C^n(z)G(z) + G(z). (3.5) 



Using (j3_4f, 

u{z) - u{z){G{z) -I) = G{z) - I. (3.6) 



Definition 3.2. We use C[G\T] to refer to the operator (bounded on L (Y) provided G € L°°(r)j in (3.6). 

In what follows we assume at a minimum that the RHP is well-posed, or C[G; T]^ 1 exists and is bounded 
ini 2 (r) and G-/G L 2 (r). 

We need to establish the smoothness of solutions of (3.6) since we approximate solutions numerically. 
This smoothness relies on the smoothness of the jump matrix in the classical sense along with a type of 
smoothness at the self-intersection points, 70, of the contour. 

Definition 3.3. Assume a £ 70- Let r 1; ...,r m be a counter-clockwise ordering of subcomponents of T 
which contain z — a as an endpomt. For G € W k '°°(T) n H k (T) we define Gi by G| f . if Y l is oriented 
outwards and (GlrJ -1 otherwise. We say G satisfies the (k — l)th-order product condition if using the 
(k — \)th-order Taylor expansion of each Gi we have 

m 

HG t =I + 0((z-a) k ),forj = l,...,k-l, Va e 7o . (3.7) 



To capture all the regularity properties of the solution of (3.6), and to aid the development of numerical 
methods, we introduce the following definition in the same vein as Definition |3.3| 



Definition 3.4. Assume that a € 70 and letTi,... ,T m be a counter-clockwise ordering of subcomponents 
ofT which contain z — a as an endpoint. For f € H k (T), define 



_ I -lim z ^ a f\ Ti (z) ifTi is oriented outward, , g ^ 

lmiz-^a (j^Y /|r i (^r) ifTi is oriented inward. 



We say that f satisfies the (k — \)th-order zero-sum condition if 

m 

fi j) = 0, for j = 0, . . . , k - 1 and Va e 7o . (3.9) 
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Figure 1: A typical intersection point and the definition of G;. 



Remark 3.2. These definitions imply conditions also when T has an isolated endpoint. In that case G = I 
at this endpoint and / = 0. 

This motivates the following definition. 

Definition 3.5. A RHP [G;T] is said to be k-regular if 

• G, G _1 satisfy the (k — 1) th- order product condition, and 

• G - 1, G- 1 - 1 e w k -°°{T) n H k (T). 

It is worth noting that when T is bounded the H k (T) condition is trivial. We include it for completeness. 
A useful result [HI Chapter II, Lemma 6.1] is: 

Lemma 3.1. Let T be a smooth non-closed curve from z — a to z — b with orientation from a to b. If 
f G H x (r), then 

D rm dC= m-m + r°m dC . 

Jr C - * a-z o-z Jv Q~ z 

Due to cancellation, Cp commutes with the weak differentiation operator for functions that satisfy the 
zcroth-order zero-sum condition. We use the notation H k (T) for the closed subspace of H k (T) consisting of 
functions which satisfy the (k — l)th order zero-sum condition. This allows us to use the boundedness of 
Cp on L 2 (r) to show the boundedness (with same norm) of Cp from H k (T) to H k (T). The results in [16] . 
combined with the fact that differentiation commutes with the Cauchy operator on these zero-sum spaces 
can be used to prove the following theorem. 

Theorem 3.1. Given a RHP [G;T] which is k-regular, assume C[G;T] is invertible on L 2 {T). Then u G 
H k (Y), the solution of (3.6), satisfies 



D k u = C[G; r]" 1 1 D (G - I) + D k {C^u ■ (G - /)) - Cp D k u ■ (G - I) ). (3.10) 



where the right-hand side of (3.10) does not depend on D u 
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Corollary 3.1. Under the hypotheses of Theorem 3.1, u satisfies an inequality of the form 



Mhht) < Pk (\\G- i-||^.~ (r) ||C[G; r]" 1 lU^cr))) \\C[G; r]" 1 |U (z , 2(r)) ||G - 7||^ (r) , (3.11) 
where pk is a polynomial of degree k whose coefficients depend on \\Cr\\c(L 2 (r))- 



Proof: Taking the norm of (3.101 gives a bound on the semi-norm ||D fe u||^2(p) in terms of ||w|| jyk-i£Q. Us- 
ing (3.101 for k — 1 gives a bound in terms of IMIjjjs-a^r). This process produces a bound of the form (3.11 ). ■ 



Remark 3.3. The expression for the derivative in Theorem \3.1\ can be used to bound Sobolev norms of the 
solution in terms of Sobolev norms of the jump matrix if a bound on the norm of the inverse operator is 
known. In some cases, when the jump matrix depends on a parameter, and the Sobolev norms of the jump 
matrix are bounded or decaying, the resulting bounds are of use. 

We often use the following theorem which is derived from the results in [7]. 

Theorem 3.2. Consider a sequence of RHPs {[G^]T]}^>q on the fixed contour T which are k-regular. 
Assume G e ^ G in W k >°°(T) n H k (T) as £ -> oo, then 

• If C[G; r] is invertible then there exists a T > such that C[G^; T] is also invertible for £ > T . 

• If is the solution of[G^,T], and $ the solution of[G;T], then $^ — $ ± -> in H k (T). 

• ||$£ — $||wi.«>(s) — > 0, for all j > 1, whenever S is bounded away from T. 

Proof: The first statement follows from the fact that C[G^;T] converges to C[G;T] in operator norm. The 
second property follows from Corollary |3.1| The final property is a consequence of the Cauchy- Schwartz 
inequality and the fact that ||ttj — w]|z,2(r) — > 0. B 



4 The Numerical Solution of Riemann— Hilbert Problems 

The goal in this section is to introduce the necessary tools to approximate the operator equation 

C[G\T]u = G-I. (4.1) 

We start with two projections I n and V n , both of finite rank. Assume V n is defined on ifj(r) and I n is 
defined on iJ 1 (r). Define X n — ranP n and Y n = ranl„ equipping both spaces with the inherited L 2 (T) 
norm. We obtain a finite-dimensional approximation of C[G;T] by defining 

C n [G;T}u=l n C[G;T}u. 

It follows that C n [G; T] : X n — > Y n . This is useful under the assumption that we can compute C[G; T] exactly 
for all u € X n . An approximate solution of (4.1 1 is obtained by 

u n =C n [G;T]- 1 X n (G-I), 

whenever the operator is invertible. We use the pair (I n ,V n ) to refer to this numerical method. 

Remark 4.1. One may have concerns about C n [G;T] because the dimension of the domain and that of the 
range are different. It turns out that C[G;T] maps H\{T) to a closed subspace o/_ff 1 (r) and we can define I n 
on this space. In the numerical framework of JTTj, solving the associated linear system results in a solution 
that must satisfy the zeroth-order zero-sum condition, justifying the theoretical construction above. In what 
follows we ignore this detail. 

To simplify notation, we define T[G; T]u = C^u(G - I) (so that C[G; T] = I - T[G; T]) and %[G; Y] = 
l n T[G;T]. We use a few definitions to describe the required properties of the projections. 



G 



Definition 4.1. The approximation C n [G;T] to C[G;T] is said to be of type (a,/3,"f) if, whenever C[G;T] is 
invertible for n > N, C n [G;T] is invertible and 

• \\C n [G;r]\\ c{Hl(r)!Yn) <C 1 n«(l+||G-7|Uoo (r) ||C I T|U( i2 (r))), 
. \\C n [G-T]-^\ c{YM <C 2 n^\\C[G ] T]-^\ c{L 2 {T)) and 

• \\Tn[G;T]\\ c(Xn ,Y n ) < C3^ 7 ||G-/|| L =e (r) ||Cp || £(L 2 (r)) . 

The constants here are allowed to depend on T. 

The first and second conditions in Definition |4.1| are necessary for the convergence of the numerical 
method. This will be made more precise below. The first and third conditions are needed to control 
operator norms as G changes. It is not surprising that the first and the third conditions are intimately 
related and in ^6] we demonstrate the connection. 

Remark 4.2. Some projections, mainly those used in Galerkin methods, can be defined directly on L 2 (T). 
In this case we replace the first condition in Definition \4-l\ with 

\\Cn[G;F]\\£,(L 2 (r),Y„) < C*in Q (l + \\G - I\\L°°(r)\\Cr\\c{L 2 (r)))- 
This condition and the second condition with a — 7 are implied by requiring 

\\Zn\\c(L 2 (r),Y n ) < C\n a . 

In this sense Galerkin methods are more natural for RHPs, though we use the collocation method of 111]/ 
below because a Galerkin method has yet to be developed. 

Definition 4.2. The pair (I n ,V n ) is said to produce an admissible numerical method if 

• The method is of type (a, /J, 7). 

• For all m > 0, \\I n u — u||#i(r) and \\VnU — u||#i(r) tend to zero faster than n~ m as n — > 00 for all 
u G H S (T); for some s. 

• I n is bounded from C(T) to L 2 (T), uniformly in n. 

Remark 4.3. We assume spectral convergence of the projections. This can be relaxed but one has to spend 
considerable effort to ensure a, (3 and 7 are sufficiently small. 

Next, we prove the generalized convergence theorem. 

Theorem 4.1. Assume that (I n ,V n ) produces an admissible numerical method. IfC[G;T] is 1-regular and 
invertible on L 2 (T), we have 

\\u - u n \\ L 2 {r) < (1 + cn a+0 )\\V n u - uII^d with (4.2) 
c = (7||C[G;r]- 1 |U (i2(r)) (l + ||G-J|U» (r) ||CF|| Ai 2 (r)) ). (4.3) 

Proof: First, for notational simplicity, define /C„ = C n [G;T], K, = C[G;T] and / = G — I . Then u n = 
K,- X l n f = K,- l l n Ku. Further, since u e Hl(T), 

u — u n =u — V n u + V n u — u n 

=u — V n u + V n u — K^InlCu 

=u - V n u + JC^JCnVnU — K~ x T n Ku 

=u - V n u + K.^ 1 (K n V n u - l n Ku) 

=U - V n U + IC^lnJClVnU - u). 

We used JC n V n u = I n K,V n u for u e Hl(T) in the last line. Taking an L 2 (T) norm, we have 

\\u- u n ||i,2(r) < \\{I + K-n 1 ' I nK-{u-'P n u)\\ L 2 { Y ) < (1 + cn a+ P)\\u ~ V n u\\ H i-(r)-* 
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Remark 4.4. In the case mentioned above where X n and V n can be defined directly on L 2 (T) we obtain a 
purely L 2 (T) based bound 

\\u - u n \\ L 2 (r) < (1 + cn a+f3 )\\u - V n u\\ L 2 {v) . 



Corollary 4.1. Under the assumptions of Theorem J^.l and assuming that [G; T] is k -regular for large k 
(large is determined by Definition 4-2) we have that $„ = I + Cru n is an approximation of the solution 
of [G; r] , in the following sense. 

• $± - $± -> in L 2 (T) and 

• ||$„ — ^H^/j.ooj^g) — > for all j > 0, whenever S is bounded away from T. 

Proof: The first claim follows from the boundedness of the Cauchy operator on L 2 (T) and, as before, the 
Cauchy-Schwartz inequality gives the second. I 

Below we always assume the numerical method considered is admissible. The ideas presented thus far 
are general. In specific cases the contour T consists of disjoint components. We take a different approach to 
solving the RHP in this case. 

Example 4.1. Consider the RHP [G; T] with T = T\ U T2 where T\ and T2 are disjoint. To solve the full 
RHP, we first solve for 3>i — the solution of [Glr^Ti] - assuming that this sub-problem has a unique 
solution. The jump on T2 is modified through conjugation by $1. Define 

G 2 = $iG|r 2 *r 1 - 

The solution $2 of ^2^2] is then found. A simple calculation shows that $ = $i$2 solves the original 
RHP [G; T] . 

This idea allows us to treat each disjoint contour separately, solving in an iterative way. When using this 
algorithm numerically, the dimension of the linear system solved at each step is a fraction of that of the full 
discretized problem. This produces significant computational savings. We now generalize these ideas. 

Consider a RHP [G; T] where T = Ti U • • • U Tg. Here each Ti is disjoint and 1^ = ctjfij + /3j for some 
contour f2j. We define Gi(z) — G(z)\r i and Hi(k) = Gj(o,fc + /3,). As a notational remark, we always 
associate Hi and G is this way. 

Remark 4.5. The motivation for introducing the representation of the contours in this fashion will be made 
clear below. Mainly, this formulation is important when at and/or /3j depend on a parameter but Qi does 
not. 

We now describe the general iterative solver. 
Algorithm 4.1. (Scaled and Shifted RH Solver) 

1. Solve the RHP [Hi;Q,i] to obtain $1. We denote the solution of the associated SIE as U\ with domain 



Qi. Define $i(z) = $1 



2. For each j = 2, . . . , £ define $j = $i(ctjZ + /3,-) and solve the RHP [Hj] Qj] with 



If, % :.,•••■!>;. •••<!" 1 



to obtain Again, the solution of the integral equation is denoted by Uj with domain Qj. Define 



3. Construct 3? = $£••• $i, which satisfies the original problem. 
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When this algorithm is implemented numerically, the jump matrix corresponding to Hj is not exact. It 
depends on the approximations of each of the <J>i for i < j and more specifically, it depends on the order of 
approximation of the RHP on fij for i < j. We use the notation rii — (ni, . . . , n*) where each rii is the order 
of approximation on fij. We use n > m whenever the vectors are of the same length and rij > rrij for all j. 
The statement n — ¥ oo means that each component of n tends to oo. Let ^ij^rii be the approximation of 
&i t j and define 

Hj, n . = <\>, i .;.*//. • • • $i,j,m-H j ^r,i,n 1 ' ' ' <h . !,../(. • 

If the method converges then Hj n . — > Hj uniformly as rij — > oo. 

A significant remaining question is: "how do we know solutions exist at each stage of this algorithm?" 
In general, this is not the case. C[G; T] can be expressed in the form JC — T where JC is the block-diagonal 
operator with blocks C[G;;Ti] and T is a compact operator. Here T represents the effect of one contour on 
another and if the operator norm of T is sufficiently small solutions exist at each iteration of Algorithm |4.1| 
This is true if the arclength of each I\ is sufficiently small. We leave a more thorough discussion of this to 
§5.1| An implicit assumption in our numerical framework is that such equations are uniquely solvable. 

The final question is one of convergence. For a single fixed contour we know that if (I n ,V n ) produces 
an admissible numerical method and the RHP is sufficiently regular, the numerical method converges. This 
means that the solution of this RHP converges uniformly, away from the contour it is defined on. This 
is the basis for proving that Algorithm |4.1| converges. Theorem |3.2| aids us when considering the infinite- 
dimensional operator for which the jump matrix is uniformly close, but we need an additional result for the 
finite-dimensional case. 

Lemma 4.1. Consider a sequence of RHP 's {[G^;T]}^>o on the fixed contourT which are k-regular. Assume 
— >• G in L°°(T) n L 2 (T) as £ — >• oo and [G; T] is k-regular, then 



If C n [G]T] is invertible, then there exists T{n) > such that C n [G^]T\ is also invertible for £ > T(n). 



• If ^n,{ is the approximate solution of [Gt; T] and is the approximate solution of [G; T], then 

— > in L 2 (T) as £ — ^ oo for fixed n. 

• ||$n,$ — &n\\wi-°°(s) 0, as £ — > oo, for all j > 1, whenever S is bounded away from T for fixed n. 
Proof: We consider the two equations 

C n [G^;T}u n ^ =I„(G 4 - I), 
C n [G;T]u n =l n {G-I). 



Since the method is of type (a, /3, 7), we have (see Definition (4.1 1), 

\\C n [G^T] - C n [G-T]\\ c(XM < C 3 nlCf |U (ia(r)) ||G 4 - G\\ L ~ {T) = E(£)rt 
For fixed n, by increasing £, we can make E(£) small, so that 

\\C n [G^r}-C n [G;T}\\ c{XnXn} < i- 1 n~ p < ' 



G 2 ||C[G;r]-i|| £(i2(r)) " \\C n [G;r}-i\\ c{ y M 
Specifically, we choose £ small enough so that 

2G 2 G 3 ||C[G;r]-i|| £(L2(r)) 



Using Lemma 2.1 C n [G^;T] is invertible, and we bound 

WCnlG^T}- 1 -^[Gjr]- 1 ^,^) < 2C 2 n 2 ^\\C[G ] T]- 1 \\l (L , (T)) E{^. (4.4) 
Importantly, the quantity on the left tends to zero as £ —¥ oo. We use a triangle inequality 

\\u„-u nti \\ L 2 (r) < ||(C n [G ? ;r]- 1 -C„[G;r]- 1 )I„(G~/)|| i 2 (r) + ||C„[G;r]- 1 l-„(G-G ? )|U 2(r) . 
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Since we have assumed that T is bounded and that the norm of I n : C(T) — > L 2 (T) is uniformly bounded in 
n, we obtain L 2 convergence of u n to u n £ as £ — > oo: 

IK " u»,*I|l»(T) < C 3 n 2 ^E(0\\G - I\\l°°(t) + C 4 n p \\G- G e || L oo (r) < C 5 n 2 ^E(0- (4.5) 
This proves the three required properties. ■ 

Remark 4.6. A good way to interpret this result is to see E(£) as the difference in norm between the 
associated infinite- dimensional operator which is proportional to the uniform difference in the jump matrices. 



Then (4.4) gives the resulting error between the finite- dimensional operators. It is worthwhile to note that if 



a — (3 = 7 = then S can be chosen independent of n. 

Now we have the tools needed to address the convergence of the solver. We introduce some notation to 
simplify matters. At stage j in the solver we solve a SIE on Clj. On this domain we need to compare two 
RHPs: 

• [Hj; rij] and 

• ll. n. 

Let Uj be the exact solution of this SIE which is obtained from [Hj]Qj]. As an intermediate step we need 
to consider the numerical solution of [Hj\Q.j\. We use Uj tn . to denote the numerical approximation of Uj 
of order rij. Also, Uj^rij will be used to denote the numerical approximation of the solution of the SIE 
associated with [-f/^rij ■; Qj]- 

Theorem 4.2. Assume that each problem in Algorithm ^. 1\ is solvable and k -regular for sufficiently large k. 
Then the algorithm converges to the true solution of the RHP. More precisely, there exists AT, such that for 
rii > Ni we have 

\\Ui, ni - Ui\\ L 2 (Qi) < C k [(maxn 4 ) Q+/3 + (maxn,) 2 ^ 7 ] ' max \\l n U, - Uj\\ H x {aj) , 
where I n is the appropriate projection for ftj . 

Proof: We prove this by induction. Since Ui_n x — Ui, ni the claim follows from Theorem |4.1| for i = 1. 
Now assume the claim is true for all j < i. We use Lemma |4. 1| to show it is true for i. Using the triangle 
inequality we have 

\\Ui,rii - E^IU 2 (n,) < \\U t . n% - t^nJU 2 ^) + \\Ui - Ui <ni \\ L 2^ Q .y 



Using Theorem 4.1 we bound the second term: 

\\Ui - C/UMnd ^ Cn^\\I n Ui - U^^y 



To bound the first term we use (4.5) 



\\Ui, nt ~ U i>ni \\ L 2 m < Cnf +1 E{rn^). (4.6) 

_E(rii_i) is proportional to the uniform difference of Hi and its approximation obta ined through the numerical 
method, H% t rii-x ■ By the induction hypothesis, if k is sufficiently large, Lemma 4. 1 tells us that this difference 



tends to zero as n^i — > oo, and the use of (4.5 ) is justified. More precisely, the Cauchy-Schwartz inequality 
for each fij, j < i and repeated triangle inequalities results in 

\\Hi - iWJU~(n,) < CJ2\\Uj - U jinj \\ L , {n . y (4.7) 



Combining (4.6) and (4.7) we complete the proof. 



Remark 4.7. The requirement that k is large can be made more precise using Definition \4-.£\ with m = 
max{Z(2a + 7), l(a + /3)}. There is little restriction if (a, f3,j) — (0,0,0). 
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5 Uniform Approximation 



In this section we describe how the above results can be used. First, we briefly describe how to obtain an 
explicit asymptotic approximation. Second, we use the same ideas to explain how numerics can be used to 
provide asymptotic approximations. The idea we continue to exploit is that the set of invertible operators 
between Banach spaces is open. Before we proceed, we define two types of uniform approximation. Let 
{^n}?>o be a sequence, depending on the parameter £, in a Banach space such that for each £, ||J7| — U^\\ — > 
as n — > oo for some f/s 

Definition 5.1. We say the sequence {t/|}f>o is weakly uniform if for every e > there exists a function 
N(£) : M + — > N taking finitely many values such that 

- Ut\\ < e. 

Definition 5.2. We say the sequence {U^}^>q is strongly uniform (or just uniform J if for every e > there 
exists N G N such that for n > N 

\\u*-ue\\<e. 

The necessity for the definition of a weakly uniform sequence is mostly a technical detail, as we do not 
see it arise in practice. To illustrate how it can arise we give an example. 

Example 5.1. Consider the sequence 

{U*}n,i>o = Isinf + e-" 2 +e-(«-"> 2 } 

I. J n,£>0 

For fixed £ ; C/| — > sin£. We want, for e > 0, while keeping n bounded 

|£/«-sin£| = | e - n2 +e-^-"> 2 | <e. 

We choose n > £ or if £ is large enough we choose < n < £. To maintain error that is uniformly less then 
e we cannot choose a fixed n; it must vary with respect to £. When relating to RHPs the switch from n > £ 
to < n < £ is related to transitioning into the asymptotic regime. 



5.1 Nonlinear Steepest Descent 

For simplicity, assume the RHP [G^,r^] depends on a single parameter £ > 0. Further, assume 
a(£)f2 + where SI is a fixed contour. It follows that the matrix-valued functions 

£r*(fc) = G*(a(Ofi + j8(0), 
are defined on Q. Assume there exists a matrix if such that \\H — if^||£»°(n) — > as £ — > oo. Use $^ and 



<& to denote the solutions of [if^;f2] and [if;f2], respectively. Applying Theorem 3.2 approximates $ in 
the limit. If everything works out, we expect to find 

^ = 9 + 0(C v ), v>{). (5.1) 

Furthermore, if we find an explicit solution to [if; f2] we can find an explicit asymptotic formula for the RHP 
[G* ; r^] . This function $ that solves [if^ ; fi] asymptotically will be referred to as a parametrix. 

We return to Example |4.1| and introduce a parameter into the problem to demonstrate a situation in 
which the RHPs on each disjoint contour decouple from each other. 

Example 5.2. Assume = rf U r| and 

rf =a x (0^1+ /3i, 

r| =a 2 (£)«2 + /32, |/?i-/3 2 |>0. 
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Assume that each fli is bounded. We consider the L 2 norm of the Cauchy operator applied to a function 
defined on and evaluated on r|: C r su(z)\ r z . Explicitly, 



1 f u(s) 



This is a Hilbert- Schmidt operator, and 



A simple change of variables shows that 



l C relli (i2(r!))L3(ri)) < / \aUe) S -J(tl + 8i-8,\*- ^ 



Since the denominator in the integral in (5.3) is bounded away from zero and both Qi are bounded, the 
right-hand side tends to zero if either ot\ or 0L2 tend to zero. 

This argument, with more contours, can be used to further justify Algorithm ^. 1\ in this limit by noting 
that this type of Cauchy operator (evaluation off the contour of integration) constitutes the operator T in fj^J 
We have the representation 



c[G« ; r«] = 



C[G\;T{] 
C[G«;T«] 



where ||T ? |U(z»(r)) ^ as ^ 00. 

This analysis follows similarly in some cases when Bi depends on £. For example, when 

M\f3i(t)-8 2 (0\=6>0. (5.4) 



One can extend this to the case where (5.4) is not bounded away from zero but approaches zero slower than 
ai(£)a 2 (£). For simplicity we just prove results for Bi being constant. 

Furthermore, the norms of the inverses are related. When each of the aj(£) are sufficiently small, there 
exists C > 1 such that 



1 

C 



\\C[G; iri^r)) ^ max rr '^(^(r,)) < C\\C[G; r]" 1 ^^^)). (5.5) 



Due to the simplicity of the scalings we allow the norms of the operators C[Gi \ I\] and C[Gi \ T^] -1 to coincide 
with their scaled counterparts C[Hi\ Oj] and C[Hi\ f2j] _1 . 

The choice of these scaling parameters is not a trivial task. We use the following rule of thumb: 

Assumption 5.1. If the jump matrix G has a factor and Bu corresponds to a qth order stationary point 
Bk i.e., 9{z) ^ C(z — Bk) q ), then the scaling which achieves asymptotic stability is ctk(Q = l^l" 1 ^ ■ 

We prove the validity of this assumption for the deformations below on a case-by-case basis. We need 
one final result to guide deformations. We start with a RHP posed on an unbounded and connected contour. 
In all cases we need to justify the truncation of this contour, hopefully turning it into disconnected contours. 

Lemma 5.1. Assume [G;T] is k-regular. For every e > there exists a function G e defined on T and a 
bounded contour T e C T such that: 

• G e = i on r \ r e , 

• \\G e — G|| 1 y<«.°°(r)nff fc (r) < e 
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• [G e ,r e ] is k -regular and 

\\C[G;T] - C[G e , r]|| £(i 2 (r)) < e||C~ \\c(L'(r))- 
Proof: A matrix-valued function / is chosen such that 

. f\ Vi e c°°(ri), 

• / = I in a neighborhood of all intersection points, 

• / has compact support, and 

• — I)f + 1 — G\\ W h,c«^ nH k(p^ < e. 

Equate G e = (G — I) f + I . The last property follows immediately. ■ 

This justifies the truncation of infinite contours to finite ones and it shows this process preserves smooth- 
ness. When it comes to numerical computations, we truncate contours when the jump matrix is, to machine 
cpsilon, the identity matrix. In what follows we assume this truncation is performed and we ignore the error 
induced. 

5.2 Direct Estimates 

As before, we are assuming we have a RHP that depends on a parameter £, [G^; T^] and is bounded. Here 
we use a priori bounds on the solution of the associated SIE which are uniform in £ to prove the uniform 
approximation. In general, when this is possible, it is the simplest way to proceed. 

Our main tool is Corollary |3.1[ We can easily estimate the regularity of the solution of each problem 
[ifj;f2j] provided we have some information about ||C[Gf ; rf ] — 1 || J c(i 2 (r i )) or equivalently \\C[Hf ■,^li]^ 1 \\r,(L 2 (fli))■ 
We address how to estimate this later in this section. First, we need a statement about how regularity is 
preserved throughout Algorithm |4.1| Specifically, we use information from the scaled jumps Hi and the local 
inverses C\Hi\ f2i] -1 to estimate global regularity. The following theorem uses this to prove strong uniformity 
of the numerical method. 

Theorem 5.1. Assume 

• {[G^,r^]}^>o is a sequence of k-regular RHPs, 

• the norm of C[Hf , f^] -1 is uniformly bounded in £, 

• \\H', [\\w k '°°(ni) < and 

• a,(£) — > as £ — > oo. 

Then if k and £ are sufficiently large 



Algorithm 4-1 applied to {[G^,r^]}^> has solutions at each stage, 



\\Uj\\H»(pi) < p h where P k depends on \\H^\\ Hk(n . )nWk ,^ (n%) , \\C[Hl;Sli] 1 \\ C (L 2 (n l )) and \\U^\\^ nj ) 
for j < i and 



• the approximation U^ n . of Iff (the solution of the SIE) at each step in Algorithm 
uniformly in £ as -> oo, that is, convergence is strongly uniform. 



4.1 



converges 



Proof: First, we note that since a«(£) — > 0, (5.3) shows that jump matrix Hf for the RHP solved at stage 



in Algorithm 4.1 tends uniformly to Hf. This implies the solvability of the RHPs at each stage in Algorithm 



4.1[ and the bound 
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for sufficiently large £. As before, C can be taken to be independent of £. We claim that ||C^||_ff»(n i ) is 
uniformly bounded. We prove this by induction. When i = 1, C7J = C[fl|,f2j] _1 (if. — J) and the claim 
follows from Corollary |3.1| Now assume the claim is true for j < i. All derivatives of the jump matrix 



Hf depend on the Cauchy transform of U§ evaluated away from fij and Hf. The former is bounded by 
the induction hypothesis and the later is boun ded by assumption. Again, sing Corollary |3.1| we obtain the 



uniform boundedness of \\llf II Theorem 



4.2 



implies that convergence is uniform in £. 



The most difficult part about verifying the hypotheses of this theorem is establishing an estimate of 
|| C[ff| ; Of ] _1 \\c(L 2 (ni)) as a function of £. A very useful fact is that once the solution of the RHP [G^; 
is known then the inverse of the operator is also known: 

C[G^}~ l u = C+ [«(* € ) + ][(* £ )~ 1 ] + ~ C- ( [u(^)+}[(^)- 1 }- ■ (5.6) 

This is verified by direct computation. When ^ is known approximately, i.e., when a parametrix is known, 
then estimates on the boundedness of the inverse can be reduced to studying the L°° norm of the parametrix. 
Then (5.5 ) can be used to relate this to each C{G\ ; rf] _1 which gives an estimate for the norm of C[Hf\ Qj] _1 . 



We study this further in section [7] 
5.3 Failure of Direct Estimates 

We study a toy RHP to motivate where the direct estimates can fail. Let <j>(x) be a smooth function 
with compact support in (— 1, 1) satisfying maxr_x i] \4>( x )\ — 1/2- Consider the following scalar RHP for a 
function /x: 

p + (x) = vT(x)(l + 0(z)(l + r 1/2 e^)), (5.7) 

/i(oo) = 1, x G [-1,1], £ > 0. (5.8) 

This problem can be solved explicitly but we study it from the linear operator perspective instead. From 
the boundedness assumption on (f>, a Neumann series argument gives the invertiblity of the singular integral 



operator and uniform boundedness of the L inverse in £. Using the estimates in Corollary 3.1 we obtain 



useless bounds, which all grow with £. Intuitively, the solution to (5.7) is close, in L to the solution to 



v + {x) =v-(x)(l + <j){x)), (5.9) 
v(oo) = 1, x e [-1,1], £> 0, (5.10) 

which trivially has uniform bounds on its Sobolev norms. In the next section we introduce the idea of a 
numerical parametrix which helps resolve this complication. 

5.4 Extension to Indirect Estimates 

In this section we assume minimal hypotheses for dependence of the sequence {[G^; r^]}^>o on £. Specifically 
we require only that the map £ i— > £/f is continuous from R + to L°°(Oj) for each i. We do not want to 
hypothesize more as that would alter the connection to the method of nonlinear steepest descent which only 
requires uniform convergence of the jump matrix. In specific cases, stronger results can be obtained by 
requiring the map £ H> fl| to be continuous from K + to W k '°°(fli). 

The fundamental result we need to prove a uniform approximation theorem is the continuity of Algorithm 



4.1 with respect to uniform perturbations in the jump matrix. With the jump matrix G we associated Hj, 
the scaled restriction of G to Tj. With G we also associated Uj, the solution of the SIE obtained from 
[Hj;Qj]. In what follows we will have another jump matrix J and analogously we use Kj to denote the 
scaled restriction of J and Pj to denote the solution of the SIE obtained from [l^,fy. 

Lemma 5.2. Assume {[G^, r^]}^>o *s a, sequence of 1-regular RHPs such that £ M- Hf is continuous from 
R + to L°°(f2j) for each i. Then for sufficiently large, but fixed, rii the map £ i— > Uf n . is continuous from 
M+ to L 2 (Oi) for each i. 
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Proof: We prove this by induction on i. For i = 1 the claim follows from Lemma 4.1 Now assume the 
claim is true for j < i and we prove it holds for i. We show the map is continuous at r/ for r\ > 0. First, from 
Lemma 14.11 



m nt -Ul n J\ L > {n) <Cn* a+ -<E(t, V ), 
where E(£,ri) is proportional to \\H^ ni i — Hf Tli _ i \\i J r X ,^f li y A similar argument as in Theorem 



4.2 



gives 



3=1 

for |£ — rj\ sufficiently small. By assumption, the right-hand side tends to zero as £ — ¥ r\ and this proves the 
lemma. ■ 



It is worthwhile noting that the arguments in Lemma |5.2| show the same continuity for the infinitc- 
dimcnsional, non-discretized problem. Now we show weak uniform convergence of the numerical scheme on 
compact sets. 

Proposition 5.1. Assume {[G^,r^]}^>o is a sequence of k -regular RHPs such that all the operators in 
are invertible for every £. Assume that k is sufficiently large so that the approximations from 
converge for every £ > 0. Then there exists a vector-valued function N(i,£) that takes finitely 



4.1 



Algorithm 
Algorithm 
many values such that 



4.1 



^AT(«)-^il'- ' 



< e. 



Moreover if the numerical method is of type (0,0,0) then convergence is strongly uniform. 



5.2 



that the function E(£, n, i) 



\U 



-u 



£ 2 (r\ 



Proof: Let S C K + be compact. It follows from Lemma . 

is a continuous function of £ for fixed n. For e > find an such that E(£,nt,i) < e/2. By continuity, 
define <5 4 (n ? ) > so that E(s,n^,i) < e for |s-f| < Define the ball B(£, S) = {s G R+ : \s-\\ < 6}. The 
open sets {-B(£, 5^)}^s cover S and we can select a finite subcover {B(£j, )}jLi- We have E(s, n^,i) < e 
whenever s € £?(£j,#e.). To prove the claim for a method of type (0,0,0), we use the fact that 5^ can be 
taken independent of and that E(s,n,i) < e for every n > n^. ■ 



Definition 5.3. Given a sequence of k-regular RHPs {[G^,r^]}^>o such that 

• r« =rju---urf, and 

another sequence of k-regular RHPs {[J^,E^]}^> is said to be a numerical parametrix if 



• = £f U • • • U 



v 



• For a^/ z 

uniformly on f2j as £ — > oo, 



J £ (7i(0* + ffi) - G«(a 4 (0fc + A) -> 0, 



(5.11) 



</ie norms of the operators and inverse operators at each step in Algorithm \4-l\ are uniformly bounded 
in £, implying uniform boundedness of in £, and 



• the approximation P^ n . of P^ (the solution of the SIE) at each step in Algorithm 
formly as minn, — > oo. 



converges uni- 
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This definition hypothesizes desirable conditions on a nearby limit problem for the sequence {[G^, r^]}^>o- 
Under the assumption of this nearby limit problem we are able to obtain a uniform approximation for the 
solution of the original RHP. 

Lemma 5.3. Assume there exists a numerical parametrix { Jj, E^}^>o for a sequence of RHPs {[G^, r^]}j>o- 
Then for every e > there exists Ni and T > such that, at each stage in Algorithm^. 1\ 



N, 



-Uf\\ L 2 m <efor^>T. 



(5.12) 



Furthermore, if the numerical method is of type (0,0,0), then (5.12) is true with Ni replaced by any Mi > 



Proof: At each stage in Algorithm 4.1 we have 

~ Etflbflio < Win* - ^njb(^) + II ^ - PfWvm + \\P$- uf\\ifl(fly (5-13) 

Since Pi n originates from a numerical parametrix we know that \\Pi n . — P^Wl 2 !^) ~^ uniformly in £ as 
rii is increased. Furthermore, \\P^ — f/ 4 - ||x,*(n) depends only on £ and tends to zero as £ — > oo. The main 
complication comes from the fact that a bound on \\Uf n . — P- 



i,Tli lli 2 (Q») 



from (4.4) depends on both rii 



and £ if the method is not of type (0,0,0). The same arguments as in Lemma 5.2 show this tends to zero. 
Therefore we choose rii large enough so that the second term in (5.13) is less than e/3. Next, we choose £ 
large enough so that the sum of the remaining terms is less than 2/3e. If the method is of type (0, 0, 0) this 
sum remains less than e if is replaced with n for n > rii. This proves the claims. I 

Now we prove the uniform approximation theorem. 

Theorem 5.2. Assume {[G^,r^]}^>o is a sequence of k-regular RHPs for k sufficiently large so that Algo- 
rithm^^ converges for each £. Assume there exists a numerical parametrix as £ — > oo. Then Algorithm ^. 1\ 
produces a weakly uniform approximation to the solution o/{[G^,r^]}^>o- Moreover, convergence is strongly 
uniform if the method is of type (0, 0, 0). 



Proof: Lemma 5.3 provides an M > and N\(i) such that if £ > M then 

lie/ 4 



,N 1 {i) 



Uf\ 



< e, for every i. 



According to Theorem 5.1 there is an N 2 (t;,i) such that 



\U' 



>r a «,o ~ ^ll i2 (°») < e ' for every l - 



The function 



N 1 (i), if£>M, 
JV a (£,i), if£<M, 



satisfies the required properties for weak uniformity. Strong uniformity follows in a similar way from Lemma 
Eland Theorem EU ■ 



Remark 5.1. This proves weak uniform convergence of the numerical method for the toy problem introduced 
in §5.5j ' we can take the RHP for v as a numerical parametrix. 

The seemingly odd restrictions for the general theorem are a consequence of poorer operator convergence 
rates when n is large. A well-conditioned numerical method does not suffer from this issue. It is worth 
noting that using direct estimates is equivalent to requiring that the RHP itself satisfies the properties of a 
numerical parametrix. 

In what follows, we want to show a given sequence of RHPs is a numerical parametrix. The reasoning for 
the following result is two-fold. First, we hypothesize only conditions which are easily checked in practice. 
Second, we want to connect the stability of numerical approximation with the use of local, model problems 
in nonlinear steepest descent. 
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Proposition 5.2. Assume 

• {[J^,E^]}j> is a sequence of k -regular RHPs, 

• C[-K"|, f2j] _1 has norm uniformly bounded in £, 

• \\ K i\\w k ^(ni) < c > and 

• 7i(£) ~~ as £ — > oo. 

TTien i/ fc and £ are sufficiently large 



Algorithm 4-1 applied to {[J^,E^]}^>o /«as solutions at each stage and 



{[J^,S^]}^>q satisfies the last two properties of a numerical parametrix (Definition 5.3) 



Proof: The proof of this is essentially the same as Theorem 5.1 



Remark 5.2. Due to the decay ofji, the invertiblity of each of C[K^,ili\ is equivalent to that of C[G^;T^]. 

This proposition states that a numerical parametrix only needs to be locally reliable; we can consider 
each shrinking contour as a separate RHP as far as the analysis is concerned. 



6 A Numerical Realization 

In [11] , a numerical framework was constructed for computing solutions to RHPs, based on a method used 
in [T2] for computing solutions to the undeformed Painleve II RHP. This framework is based on Chebyshev 
interpolants. Consider the RHP [G;T], T = Ti U ■ • • U 1^ , where each I\ is bounded and is a Mobius 
transformation of the unit interval: 

Afi([-1, 1]) = IV 

Definition 6.1. The Chebyshev points of the second kind are 



x 



[-1,1]," _ 



-1,1], n 



-1,1], n 



-1 



COS 7T ( 1 



cos ■ 



The mapped Chebyshev points are denoted 

x l > n = Afitsel- 1 ' 1 ]). 

Given a continuous function fa defined on Ij we can find a unique interpolant at x l ' n using mapped 
Chebyshev polynomials. Given a function, / defined on the whole of T, we define T n to be this interpolation 
projection applied to the restriction of / on each Ti. Clearly, 

and, because x^ n contains all junction points, 

l n :Hl{T)^Hl{T). 

The framework in |11] is given by the pair (I„,I„) and the matrix C n [G; T] is equal to T n C[G\ T]I n with 
some unbounded components subtracted; obeying the requirement that the two operators agree on ifj(r). 
Now we address the properties required in Definition |4.2| 
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Lemma 6.1. When G E VF 1,00 (r), the numerical method in fllf satisfies: 

• I n is uniformly bounded in n from C(T) to L 2 (T) when T is bounded. 

• ||C„[G;r]|| £(// i (r )^ 7i) < C(l + ||G - 7|| L oo (r) ||C r 7||£( i 2( r)) ). 

• \\r n [G;T]\\ c{ x n ,Y n ) <Cn 2 ||G-7|| L » (r) ||Crl| £ (L2 (r)) . 

• \\l n u- u\\ H i( r ) <C s n 2 - s \\u\\ H »(r)- 

Proof: First, note that these constants depend on T. Using the Dirichlet kernel one proves that T n is 
uniformly bounded from G(T) to an L 2 space with the Chebyshev weight P3- The norm on this weighted 
space dominates the usual L 2 {Y) norm, proving the first result. For the second statement we take u E Hl(T) 
and consider 

\\I n -I n {G - I)Cp M|| L 2 (r) < ||2"„||£(c(r),L 2 (r))(l + j|G - /||x,°=(r) ||Cp |U(fl-i(r),jJ 1 (r))||'"||H 1 (r))- 

Since Y n is equipped with the L 2 {T) norm and \\C^ WciHKr)^ 1 ^)) = \\Cr\\c(L 2 (r)) we obtain the second 
property. We then use for u E X n (see <j4j) that ||u||iji(r) < Gn 2 1 1 it || £2 (p) to obtain 

\\T n [G;T]\\ c(XntYn) <Gn 2 ||G-/|| L ^ (r) ||Cr||£ (L 2 (r)) . 

The last statement follows from estimates in [TJ] for the pseudo-spectral derivative. ■ 

The final property we need to obtain an admissible numerical method, the boundedness of the inverse, is 
a very difficult problem. We can verify, a posteriori, that the norm of the inverse does not grow too much. 
In general, for this method, we see at most logarithmic growth. We make the following assumption. 

Assumption 6.1. For the framework in fllf we assume that whenever [G;T] is 1-regular and C[G;r] _1 
exists on L 2 {T) as a bounded operator we have for n > N 



lic^r]- 1 !!^^) < c^||c[G ; r]- 1 || £(i2(r)) , p > o. 



(6.1) 



Remark 6.1. This assumption is known to be true for a similar collocation method on the unit circle using 
Laurent monomials [13]. 

With this assumption the numerical method associated with (I n ,l n ) is of type (0,/3,2). In light of 
Theorem |4.1 1 we expect spectral convergence and the bound in Assumption |6.1| does not prevent convergence. 
We combine Assumption |6.1[ Theorem |4. 1| and Theorem |3.1| to obtain 

\\u-u n \\ L 2 {r) < C(||C[G;r]- 1 ||£ :(i2(r)) (l + || G - (r) ||£: (i2(r)) )« 2 +^- fc ||w|| ffMr) . (6.2) 



7 Application to Painleve II 



We present the RHP for the solution of the Painleve II ODE ( |1.1[ ). Let T = Ti U • • ■ U Fk with r. ( = 
| se i7r(j/3-i/6) . fi g E+}, i.e., r consists of six rays emanating from the origin, see Figure pi The jump 
matrix is defined by G(A) = G;(A) for A E Ti, where 



Gi{x; A) = Gi(A) 



1 o --i8/3A -2teA 

o 1 

1 ' 

s . e i8/3A 3 +2ia;A ^ 



if i is even, 
if i is odd. 



From the solution $ of [G; T], the Painleve function is recovered by the formula 

u(x) = lim z<& 12 (z), 
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- ^z' s -2ixz 



1 

-ifz 3 ~ 2ixz -i 



1 

Figure 2: The contour and jump matrix for the Painleve II RHP. 



where the subscripts denote the (1, 2) entry. This RHP was solved numerically in [T2"] . 

For large |x|, the jump matrices G are increasingly oscillatory. We combat this issue by deforming the 
contour so that these oscillations turn to exponential decay. To simplify this procedure, and to start to mold 
the RHP into the abstract form in fjij wc first rescale the RHP. If wc let z = y^\x\X, then the jump contour 
r is unchanged, and 



$+(z) = $+(VMA) = $-(y/\x\\)G(y/\x\\) = $-(z)G(z), 
where G(z) — Gi(z) on I\ with 



Gi(z) 



1 s,e 



-&{*) 



£ = |x| 3 / 2 and 



1 
1 

*e* fl M 1 



2? 



if i is even, 
if i is odd, 



Then 



z) = - (4z 3 + 2e tavsx z) 



u(x) = lim A$i2(x;A) = \Jx lim z§yi(x;z). 

A— foo A— ►oo 



(7.1) 



We assume that s\s^ > 1 and x < 0. We deform T to pass through the stationary points ±1/2, resulting in 
the RHP on the left of Figure g) 
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The function 



GqG\G2 — 



1 - sis 3 



1 + S1S2 



has terms with exp(±£6>). It cannot decay to the identity when deformed in the complex plane. We can 
resolve this issue by using lensing [4j. Suppose we factor the jump function as ABC. It is possible to separate 
the single contour into three contours, as in Figure [3j assuming A and C are analytic between the original 
contour and continuous up to the new contour. If $ satisfies the jumps on the split contour, it is clear that 
we can recover 4> by defining <I> = <J>C between the top contour and the original contour, <I> — between 
the original contour and the bottom contour, and $ = $ everywhere else, see the bottom of Figure [3} As in 
the standard deformation of contours, the limit at infinity is unchanged. 





Figure 3: Demonstration of the lensing process. 



Now consider the LDU factorization: 



GqGxG 2 = LDU = 



-ce gi 



1 - S1S3 




1-S1S3 
1 



U decays near ioo, L decays near — £00, both go to to the identity matrix at infinity and D is constant. 
Moreover, the oscillators on L and U are precisely those of the original G matrices. Therefore, we reuse the 
path of steepest descent, and obtain the deformation on the right of Figure [4] The LDU decomposition is 




Figure 4: Left: Initial deformation along the paths of steepest descent. Right: The deformed contour after 
lensing. 



valid under the assumption S1S2 > 1. 
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7.1 Removing the connected contour 

Although the jump matrix D is non-oscillatory (in fact, constant), it is still incompatible with the theory 
presented in Sj5] we need the jump matrix to approach the identity matrix away from the stationary points. 
Therefore, it is necessary to remove this connecting contour. Since D = diag(c£i, (I2) is diagonal, we can solve 
P+ = P~D with P(oo) = I on (-1/2, 1/2) in closed form [7J: 



)ilogd2/27r 



PUP- 1 PG 2 U- 1 P- L 




PGsLP' 1 PL-^P- 1 

Figure 5: Left: Definition of <1> in terms of <J/. Right: Jump contour for vp. 

This parametrix solves the desired RHP for any choice of branch of the logarithm. However, we must 
choose the branch so that the singularity is square integrable [7 ■ In this case this is accomplished by choosing 
the standard choice of branch. 

We write 

$ = *P. 

Since P satisfies the required jump on (—1/2, 1/2), $ has no jump there. Moreover, on each of the remaining 
curves we have 

*+ = &+P- 1 = $"GP _1 = * _1 PGP -1 , 

and our jump matrix becomes PGP -1 . Unfortunately, we have introduced singularities at ±1/2 and the 
theory of <|5] requires some smoothness of the jump matrix. This motivates alternate definitions for \t in 
circles around the stationary points. In particular, we define $ in terms of ^ by the left panel of Figure 
[5j where has the jump matrix defined in the right. A quick check demonstrates that this definition of $ 
indeed satisfies the required jump relations. 
We are ready to apply Algorithm |4.1| Define 

fl = {z : \\z\\ = l}U{re' i7r / 4 : r G (1, 2)} U {re 3 " /4 : r € (1, 2)}U{re~ 3l7r/4 : r £ (1, 2)} U {re^/ 4 : r G (1,2)}. 

In accordance with Assumption |5 . 1 1 we have 

ri = i + r 1/2 « and r 2 = -i + r 1/2 a 

with the jump matrices defined according to Figure [5j Paths of steepest descent are now local paths of 
steepest descent. 



P{z) = 



2x+l 
2x~l 



i log di/27r 
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7.2 Uniform Approximation 

We have isolated the RHP near the stationary points, and constructed a numerical algorithm to solve the 
deformed RHP. We show that this numerical algorithm approximates the true solution to the RHP. In order 
to analyze the error, we introduce the local model problem for this RHP following [7J. 
Define the Wronskian matrix of parabolic cylinder functions £)„(£), 



Zo(0 = 



D- v -i(iQ 



d =D v {Q 



(7.2) 



and the constant matrices 

H k+2 = e i < v+1 ^ a3 H k e i ^ v+1 ^ a3 ,H 



1 " 


,H 2 = 


' 1 


hi ' 




' 1 


. h o 1 . 





1 


,0-3 = 


-1 



h 



2tt 



hi 



i> + i)' - r(- v y 

The sectionally holomorphic function Z{C) is defined as 



This is used to construct the local solutions 

* r (z) = B(z)(~h 1 /s 3 )-' 73 / 2 e lt ' J3/2 2-' T ^ 2 





if 


argCe (-tt/4,0), 


Zi«), 


if 


argCe (0,tt/2), 


< z 2 ((), 


if 


argC £ (7r/2,7r), 


Zs(C), 


if 


argCG (7r,37r/2), 




if 


argCG (37T/2, 77T/4) 



C(*) i 
i o 



Z(C(z))(~/ ll / S3 ) CT3/2 . 



tf*(z) = cr2* r (-z)cr 2 , 



where 



a 2 







1/2 
1/2 



,C(«)=2V-*fl(«)+*fl(l/2). 



Consider the sectionally holomorphic matrix- valued function "J'(z) defined by 

r p(z), if |2±i/2| > ij, 

= < tf r (z), if \z-l/2\ <R, 
{ ¥{z) 7 if \z + l/2\ < R. 

We use [G; f ] to denote the RHP solved by See the top panel of Figure[(3]for f . In [7 it is shown that W 
satisfies the RHP for $ exactly near z = 1/2 and for if? 1 near z = —1/2. Notice that ^ r and are bounded 
near z = ±1/2. In the special case where logdi el, P remains bounded at ±1/2. Following the analysis in 
[7] we write 

<f>(z) = X {z)i>{z), 

where x — > / as £ — > oo. 

Wc deform the RHP for ^ to open up a small circle of radius r near the origin as in Figure [H] We use 
[Gi;f i] to denote this deformed RHP and solution ^i. See Figure[(3]for f\. Also it follows that ^(z)P~ 1 (z) 
is uniformly bounded in z and £. Further, \&i has the same properties. Since is uniformly bounded in 
both z and £ we use (5.6 1 to show that CIGi.Yi]" 1 has uniformly bounded norm. We wish to use this to 
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show the uniform boundedness of the inverse C[G; V] . To do this we extend the jump contours and jump 
matrices in the following way. Set T e = T U Ti and define 



G e (z) = 
G e (z) = 



G{z) if z g r, 

/ otherwise, 

Gi(z) if ze f x , 
/ otherwise. 



The estimates in [7] show that G e — G e — > uniformly on T e . It follows that C[G e ;T e ] 1 is uniformly 



bounded since the extended operator is the identity operator on T\Tj_. Theorem 3.2 implies that C[G e ; r e ] _1 is 
uniformly bounded for sufficiently large £, which implies that C[G; is uniformly bounded for £ sufficiently 
large, noting that the extended operator is the identity operator on the added contours. We now use this 
construction to prove the uniform convergence of the numerical method using both direct and indirect 
estimates. 



7.3 Application of Direct Estimates 

We proceed to show that the RHP for "J satisfies the properties of a numerical parametrix. This requires 
that the jump matrices have uniformly bounded Sobolev norms. The only singularities in the jump matrices 
is of the form 

After transforming to a local coordinate k, z — £ _1 / 2 fc — 1/2, we see that 

T 1/2 fc + i 



S(k) = s(C 1/2 k - 1/2) = C lv/2 f - 



k 







8i / 










-I) 




"2 H 


-8 ZH 





The function S(k) is smooth and has uniformly bounded derivatives provided k is bounded away from k = 0. 
The deformations applied thus far guarantee that k will be bounded away from 0. To control behavior of 
the solution for large k we look at the exponent which appears in the jump matrix 

*>-!- 

and define 

e(fc) = 0(r 1/2 k - 1/2) = | - 4i fc 2 /e + 1 k 3 /e /2 - 

If we assume that the contours are deformed along the local paths of steepest descent, all derivatives of 
e £©(fc) are exponentially decaying, uniformly in £. After applying the same procedure at z = 1/2 and after 



contour truncation, Theorem |5.1| implies the RHP for ^ satisfies the hypotheses of Theorem 5T proving 
strong uniform convergence. 



7.4 Application of Indirect Estimates 

The second approach is to use the solution of the model problem to construct an numerical parametrix. Since 
we have already established strong uniform convergence we proceed to establish a theoretical link with the 
method of nonlinear steepest descent, demonstrating that the success of nonlinear steepest descent implies 
the success of the numerical method, even though the numerical method does not depend on the details of 
the nonlinear steepest descent. We start with the RHP [Gi;fi] and its solution As before, see Figure ^ 
for r\. Define u = (Wi) + — (^l) - which is the solution of the associated SIE on T±. The issue here is that 
we cannot scale the deformed RHP in Figure [5] so that it is posed on the same contour as [G; T] . We need 
to remove the larger circle. 
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Figure 6: Top: Jump contours for the model problem with solution ty. Note that J r and J/ are the jumps 
on the outside of the circles. They tend uniformly to the identity as £ —¥ oo 7\. Center: The jump contours, 
Fi, for the function The inner circle has radius r and the outer circle has radius R. Bottom: Contour 
on which U is non-zero. This can be matched up with the right contour in Figure [5J 

We define a new function U = ucf) where 4> is a C°° function with support in (B(— 1/2, R)UB(l/2, R))nTx 
such that 6 = 1 on (23(1/2, r) U B(— 1/2, r)) n Fi for r < R. Let P2 be the support of U (see bottom contour 
in Figure 6b . Define ^2=2 + Cp 2 U. From the estimates in [7] , it follows that 

= * 2 -G 2 

where G2 — G tends uniformly to zero as £ — > 00. We have to establish the required smoothness of U. We 
do this explicitly from the above expression for ^>P~ X after using the scalings z = £ -1 ' 2 fc ± 1/2. The final 
step is to let £ be large enough so that we can truncate both [G;T] and [G^;!^] to the same contour. We 



use Proposition 5.2 to prove that this produces a numerical parametrix. Additionally, this shows how the 



local solution of RHPs can be tied to stable numerical computations of solutions. 

Remark 7.1. This analysis relies heavily on the boundedness of P. These arguments would fail if we were 
to let P have unbounded singularities. In this case one approach would be to solve the RHP for \ - The jump 
for this RHP tends to the identity. To prove weak uniformity for this problem one only needs to consider the 
trivial RHP with the jump being the identity matrix as a numerical parametrix. 

7.5 Numerical Results 

In Figure [7] we plot the solution to Painleve II with (si,S2,S3) = 2, 3) and demonstrate numerically 
that the computation remains accurate in the asymptotic regime. We use u(n, x) to denote the approximate 



solution obtained with n collocation points per contour. Since we are using (7.1) we consider the estimated 



relative error by dividing the absolute error by yfx. We see that we retain relative error as x becomes large. 

Remark 7.2. Solutions to Painleve II often have poles on the real line, which correspond to the RHPs not 
having a solution. In other words, ||C[r, is not uniformly bounded, which means that the theory of 

this paper does not apply. However, the theorems can be adapted to the situation where x is restricted to a 
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subdomain of the real line such that \\C[T, is uniformly bounded. This demonstrates asymptotic stability 

of the numerical method for solutions with poles, provided that x is bounded away from the poles, similar to 
the restriction of the asymptotic formulas in [?|/. 
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Figure 7: Left: Plot of the solution, u, for small x (Solid: real part, Dashed: imaginary part). Right: Relative 
error. Solid: |x|" 1 / 2 |u(12, x) ~ u(36, x)|, Dashed: \u(8,x) -u(36,x)\/-y/\x\, Dotted: |u(4, x) - u(36, aOI/vW- 
This plot demonstrates both uniform approximation and spectral convergence. 



8 Application to the modified Korteweg— de Vries Equation 



In this section we consider the numerical solution of the modified Korteweg-de Vries equation ( |1.2[ ) (mKdV) 
for x < 0. The RHP for mKdV is [H] 

$+(z) = $~(z)G(z), zel, 
$(oo) = /, 

1 - p(z)p(-z) -p{-z)e~ e ^ 
p(z)e e ^ 1 

(z) = 2izx + 8iz 3 t. 



G(z) 



In the cases we consider p is analytic in a strip that contains R. If x <C — ci 1 / 3 the deformation is similar 
to the case considered above for Painleve II and asymptotic stability follows by the same arguments. We 
assume x = — 12c 2 i 1 / 3 for some positive constant c. This deformation is found in [15] . We rewrite 9: 

9(z) = -24ic 2 (zt 1/3 ) + 8z(zi 1/3 ) 3 . 

We note that 9 / (zq) = for zq = ±^—x/ (12i) = ±ct~ 1 / 3 . We introduce a new variable k = zt 1 ! 3 jc so that 

9(kct- 1/3 ) = -24ic 3 k + 8ic 3 fc 3 = 8ic 3 (k 3 - 3fe). 

For a function of f(z) we use the scaling f(k) = f(kct~ 1 ^ 3 ). The functions 9, G and p are identified similarity. 
After deformation and scaling, we obtain the following RHP for $(fc): 

$+0) = $-(fc) J(fc), fc g e = [-1, i] u ri u r 2 u r 3 u r 4 , 



j(k) 



1 

p{k)e 8 ^ 1 
1 -p(-k)e- S W 
1 



if fee [-1,1], 
if k e ri ur 2 , 

if A; g r 3 u r 4 , 
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where Ti, i = 1,2,3,4, shown in Figure [8j are locally deformed along the path of steepest descent. To 
reconstruct the solution to mKdV we use the formula 

u{x,t) = 2iz Q lim k$ 12 (k). (8.1) 

fc— >oo 

Remark 8.1. We assume p decays rapidly at oo and is analtytic in a strip that contains the real line. This 
allows us to perform the initial deformation which requires modification of the contours at oo. As t increases, 
the analyticity requirements on p are reduced; the width of the strip can be taken to be smaller if needed. We 
only require that each Ti lies in the domain of analticity for p. More specifically, we assume t is large enough 
so that when we truncate the contours for numerical purposes, they lie within the strip of analyticity for p. 




Figure 8: Jump contours for the RHP for mKdV. 



The parametrix derived in [5] is used to show that C [J, E] has an inverse that is uniformly bounded by 



using (5.6) as was done in the previous section. We use the analyticity and decay of p at oo along with the 
fact that the contours pass along the paths of steepest descent. 

The contour is fixed (i.e., independent of x, t and c), and this situation is more straightforward to analyze 
than the previous example. Repeated differentiation of J(k) proves that this deformation yields a uniform 
numerical approximation. Furthermore, replacing c by any smaller value yields the same conclusion. This 
proves the uniform approximation of mKdV in the Painleve region 

{(x,t):t>e, x < -e,x > -12c 2 t 1/3 }, e > 0. 

where e is determined by the analyticity of p. 



8.1 Numerical Results 

In Figurejijjwe show the solution with initial data u(x, 0) = — 2e~ x with c = \/9/4. The reflection coefficient 
is obtained using the method described in |15) . We use the notation u(n,x,t) to denote the approximate 
solution obtained with n collocation points per contour. We see that the absolute error tends to zero rapidly. 
More importantly, the relative error remains small. We approximate the solution uniformly on the fixed, 
scaled contour. When we compute the solution using (8.1 1 we multiply by zq which is decaying to zero along 
this trajectory. This is how the method maintains accuracy even when comparing relative error. 
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Figure 9: Left: Plot of the solution along x — —(St) 1 / 3 for small time. Center: Absolute error, \u(5, x, t) 
u(10,x,t)\, for long time Right: Relative error \u(5,x,t) — u(10, x, t)|/|u(10, x, t) for long time. 
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